package edu.usc.epigenome.dnaase.script;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.List;

import net.sf.samtools.SAMFileReader;
import net.sf.samtools.SAMRecord;
import net.sf.samtools.util.CloseableIterator;

import org.kohsuke.args4j.Argument;
import org.kohsuke.args4j.CmdLineException;
import org.kohsuke.args4j.CmdLineParser;
import org.kohsuke.args4j.Option;

import edu.usc.epigenome.genomeLibs.PicardUtils;




	public class SamToFpFnFilterByMinMapQ {
		/**
		 * @param args
		 */
			
			//final private static String prefix = "methylCGsRich_ASM_AllSnp_";
			final private static String USAGE = "SamToFpFnCytosine [opts] sampleName snpFileName bamFilePath ...";

			final private static int PURGE_INTERVAL = 20000; // We purge our stored Cpgs once we get this many bases past them.
			//final private static int ALLELE_GA_NUMBER = 1; //for each reads, when there are more than 1 GA position is different from reference sequence, we define it belongs to another allele. 
			
			/**
			 * @param args
			 */
			
			@Option(name="-minMapQ",usage="minimum mapping quality (default 30)")
			protected int minMapQ = 15;
			@Option(name="-minReadCov",usage="minimum read coverage (default 4)")
			protected int minReadCov = 1;
			@Option(name="-minBaseQual",usage="minimum Base quality (default 10)")
			protected int minBaseQual = 0;
			
			@Option(name="-debug",usage=" Debugging statements (default false)")
			protected boolean debug = false;
			// receives other command line parameters than options
			@Argument
			private List<String> stringArgs = new ArrayList<String>();

			
			
			/**
			 * @param args
			 */
			public static void main(String[] args)
			throws Exception
			{
				SamToFpFnFilterByMinMapQ test = new SamToFpFnFilterByMinMapQ();
				test.doMain(args);
			}

			public void doMain(String[] args)
			throws Exception {

				CmdLineParser parser = new CmdLineParser(this);
				// if you have a wider console, you could increase the value;
				// here 80 is also the default
				parser.setUsageWidth(80);
				try
				{
					parser.parseArgument(args);
					if (stringArgs.size() < 2) throw new CmdLineException(USAGE);
				
				}
				catch (CmdLineException e)
				{
					System.err.println(e.getMessage());
					// print the list of available options
					parser.printUsage(System.err);
					System.err.println();
					return;
				}

				String sampleName = stringArgs.remove(0);
				String snpFileName = stringArgs.remove(0);
				String bamFilePath = stringArgs.remove(0);
				BufferedReader br = new BufferedReader(new FileReader(snpFileName));
				String line;
				
				String preChr = "chr1";
				String fn = sampleName + "Merge_chr1.NODUPS.sorted.calmd.NODUPS.bam";
				File inputSamOrBamFile = new File(bamFilePath,fn);
				SAMFileReader inputSam = new SAMFileReader(inputSamOrBamFile);
				inputSam.setValidationStringency(SAMFileReader.ValidationStringency.SILENT);
				//String fnAllSnp = sampleName + "_SNP.filterCNV.filterSeqDepth." + minMapQ + "-" + minBaseQual + "-" + minReadCov + "-" + minAlleleCount + "-" + minAlleleFreq + "." + "txt";
				//String fnFalseSnp = sampleName + "_SNP.filterCNV.filterSeqDepth.falseSnp." + minMapQ + "-" + minBaseQual + "-" + minReadCov + "-" + minAlleleCount + "-" + minAlleleFreq + "." + "txt";
				//PrintWriter allSnpWriter = new PrintWriter(new File(fnAllSnp));
				//PrintWriter falseSnpWriter = new PrintWriter(new File(fnFalseSnp));
				String cpgSnpFileName = snpFileName + ".filterByPoorMapQ.cpg.txt";
				PrintWriter cpgWriter = new PrintWriter(new File(cpgSnpFileName));
				String cphSnpFileName = snpFileName + ".filterByPoorMapQ.cph.txt";
				PrintWriter cphWriter = new PrintWriter(new File(cphSnpFileName));
				String allSnpFileName = snpFileName + ".filterByPoorMapQ.all.txt";
				PrintWriter allWriter = new PrintWriter(new File(allSnpFileName));
				
				while( (line = br.readLine()) != null){
					String[] tmpArray = line.split("\t");
					
						
					String chr = "chr" + tmpArray[1];
					
					if(!chr.equalsIgnoreCase(preChr)){
						inputSam.close();
						fn = sampleName + "Merge_" + chr + ".NODUPS.sorted.calmd.NODUPS.bam";
						File inputNewSamOrBamFile = new File(bamFilePath,fn);
						//System.err.println(inputNewSamOrBamFile.getName());
						inputSam = new SAMFileReader(inputNewSamOrBamFile);
						inputSam.setValidationStringency(SAMFileReader.ValidationStringency.SILENT);
						preChr = chr;
					}
					
					Integer snpPosition =(int) Math.round(Double.parseDouble(tmpArray[2]));
					
					boolean cpgFlag = false;
					boolean cytosineFlag = true;
					//boolean agNegFlag = false;
					// this is just for cheat, to see the sepcific sequence depth filter to AG/TC SNP would make this kind of SNP to normal number...(normal number is 4X than other SNP) 
					//if((tmpArray[5].equalsIgnoreCase("A") && tmpArray[6].equalsIgnoreCase("G")) || (tmpArray[5].equalsIgnoreCase("T") && tmpArray[6].equalsIgnoreCase("C"))){
					//	agFlag = true;
						//System.err.println(line);
					//}
					
					if((tmpArray[5].equalsIgnoreCase("A") && tmpArray[6].equalsIgnoreCase("A")) || (tmpArray[5].equalsIgnoreCase("T") && tmpArray[6].equalsIgnoreCase("T")) ||  (tmpArray[5].equalsIgnoreCase("A") && tmpArray[6].equalsIgnoreCase("T")) || (tmpArray[5].equalsIgnoreCase("T") && tmpArray[6].equalsIgnoreCase("A")))
						cytosineFlag = false;
					int totalNum = 0;
					
					
					CloseableIterator<SAMRecord> chrIt = inputSam.queryOverlapping(chr, snpPosition, snpPosition);
					record: while (chrIt.hasNext())
					{	
						SAMRecord samRecord = chrIt.next();
						int mapQual = samRecord.getMappingQuality();
						byte[] baseQual = samRecord.getBaseQualities();
						boolean unmapped = samRecord.getReadUnmappedFlag();
						if (unmapped || (mapQual <= minMapQ))
						{
							continue record;
						}
						String seq = PicardUtils.getReadString(samRecord, true);
						try{
							
							String ref = PicardUtils.refStr(samRecord, true);
							boolean negStrand = samRecord.getReadNegativeStrandFlag();
							int alignmentS = samRecord.getAlignmentStart();
							//int readsStart = (negStrand) ? 0 - samRecord.getUnclippedEnd() : alignmentS;
							//int readsEnd = (negStrand) ? 0 - alignmentS : samRecord.getUnclippedEnd();
							int	onRefCoord = (negStrand) ? samRecord.getUnclippedEnd() : alignmentS; 
							int i = Math.abs(snpPosition - onRefCoord);
							if (seq.length() != ref.length())
							{
								System.err.println("SeqLen(" + seq.length() + ") != RefLen(" + ref.length() + ")");
								System.err.println(seq + "\n" + ref);
							}
							int seqLen = Math.min(seq.length(), ref.length());
							byte baseQS = (negStrand) ? baseQual[seqLen-1-i] : baseQual[i];
							if( baseQS <= minBaseQual )
							{
								continue record;
							}
							if (ref.charAt(i) == '-')
							{
								break record;
							}
							
							totalNum++;
							
							if(cytosineFlag){
								if(tmpArray[5].equalsIgnoreCase("C") || tmpArray[6].equalsIgnoreCase("C")){
									if((ref.charAt(i) == 'C') || (tmpArray[5].equalsIgnoreCase("C") && tmpArray[6].charAt(0)==ref.charAt(i)) || (tmpArray[6].equalsIgnoreCase("C") && tmpArray[5].charAt(0)==ref.charAt(i))){
										if(negStrand){
											if(i > 0){
												if(ref.charAt(i-1) == 'G')
													cpgFlag = true;
											}
										}
										else{
											if(i < seqLen-1){
												if(ref.charAt(i+1) == 'G')
													cpgFlag = true;
											}
										}
									}
								}
								else{
									if((ref.charAt(i) == 'G') || (tmpArray[5].equalsIgnoreCase("G") && tmpArray[6].charAt(0)==ref.charAt(i)) || (tmpArray[6].equalsIgnoreCase("G") && tmpArray[5].charAt(0)==ref.charAt(i))){
										if(negStrand){
											if(i < seqLen-1){
												if(ref.charAt(i+1) == 'C')
													cpgFlag = true;
											}
										}
										else{
											if(i > 0){
												if(ref.charAt(i-1) == 'C')
													cpgFlag = true;
											}
										}
									}
								}
							}
							
						}
						catch (Exception e)
						{
							System.err.println("-----------------------------------------");
							System.err.println("Couldn't handle seq #");
							System.err.println(seq);
							e.printStackTrace(System.err);
							System.err.println("-----------------------------------------");
//							chrIt.close();
//							System.exit(1);
						}
						
					}
					chrIt.close();
					
					
						
					if(totalNum >= minReadCov ){
						if(cpgFlag){
							cpgWriter.println(line);
						}
						else if(cytosineFlag){
							cphWriter.println(line);
						}
						allWriter.println(line);
					}

						
				}
				inputSam.close();
				cpgWriter.close();
				cphWriter.close();
				allWriter.close();
				//allSnpWriter.close();
				//falseSnpWriter.close();
				
				
				System.out.printf("over");
		}

	}